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Abstract 

A hybrid Car-Parrinello QM/MM molecular dynamics simulation has been carried 
out for the Watson-Crick base pair of 9-ethyl-8-phenyladenine and 1-cyclohexyluracil 
in deutero chloroform solution at room temperature. The resulting trajectory is an- 
alyzed putting emphasis on the N-H- ■ -N Hydrogen bond geometry. Using an em- 
pirical correlation between the N ■ • • N-distance and the fundamental NH-stretching 
frequency, the time-dependence of this energy gap along the trajectory is obtained. 
From the gap-correlation function we determine the infrared absorption spectrum 
using lineshape theory in combination with a multimode oscillator model. The ob- 
tained average transition frequency and the width of the spectrum is in reasonable 
agreement with recent experimental data. 
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1 Introduction 



The dynamics of the hydrogen bonds (HBs) linking complementary nucleic 
acid bases is an essential element responsible for the properties of DNA [1J. 
Among the various experimental methods used to study HBs, infrared (IR) 
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spectroscopy serves as an important means because of its remarkable sensi- 
tivity to HB formation and geometry [2j. In particular, absorption lineshapes 
in the HB stretching region contain information on the strength of the HB, 
its anharmonic couplings to fingerprint modes via Fermi-type resonances and 
to low-frequency intra- and intermolecular HB modes. At the same time the 
widths of the absorption profiles give evidence for the interaction with the 
surrounding medium [3] . The details of HB lineshapes as well as the responsi- 
ble dynamical processes can be investigated using nonlinear time-resolved IR 
spectroscopy [4J. So far there are only a few ultrafast IR spectroscopic studies 
of DNA samples. Zanni et al. [5] used two-dimensional IR spectroscopy in the 
carbonyl-stretching region to unravel coupling patterns sensitive to secondary 
structures in guanine-cytosine bases. A related theoretical investigation based 
on a vibrational exciton model was given by Cho and coworkers (see, e.g., Ref. 
[5]). More recently, Heyne et al. employed the two-color pump-probe method 
to identify the symmetric NH 2 stretching fundamental transition in adenine- 
thymine oligomers [7J. 

The complexity of DNA oligomers can be considerably reduced by focussing on 
the base pair building blocks. There are a number of high-resolution IR+UV 
spectroscopic studies in the gas phase by the groups of de Vries and Klein- 
ermanns. In Ref. [8], for instance, it was shown that Watson-Crick (WC) 
pairing is not the most favorable case in the gas phase (see also simulation 
in Ref. [9]). WC pairing can be enforced by substitutions as shown, e.g., for 
adenine and uracil derivatives in deuterochloroform solution (9-ethyladenine, 
1-cyclohexyluracil) [TU]. Recently, the first subpicosecond time-resolved study 
of the vibrational relaxation of the H-bonded NH stretching fundamental tran- 
sition in a substituted adenine-uracil pair has been reported by Wouterson and 
Cristalli [llj. In particular they observed that base-pairing decreases the re- 
laxation time by a factor of three as compared to monomeric uracil. Since 
vibrational energy flow in H-bonded nucleic acid bases is of apparent impor- 
tance for the understanding of DNA dynamics, not at least in the context of 
DNA photoprotection, there is a need to develop a microscopic picture for this 
process. 

In this contribution we will present the first step towards this goal which 
is the simulation of the NH stretching IR absorption lineshape around 3200 
cm -1 . In doing so we will follow the argument of Baerends and co-workers, 
that the electrostatic description of the HB in WC base pairs is questionable 
and charge transfer as well as resonance assistance occur which necessitates 
the quantum treatment of the electrons [T2]. For this purpose we employ Car- 
Parrinello Molecular Dynamics (CPMD) [T3] which has been used previously 
to simulate IR spectra, e.g., of HBs formed between water molecules and uracil 
|14j . of protonated water networks in bacteriorhodopsin [15], of intramolecular 
H-bonds in a Mannich base in gas and crystalline phase [T6|17|ll8j as well as of 
crystalline HsO^ClOj [19] and picolinic acid N-oxide [20]. In order to model 
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the effect of the fluctuating deuterochloroform environment on the N — H • • • N 
HB dynamics we have chosen to use a QM/MM hybrid method [21] as de- 
tailed in Section [2J This will provide us with a classical trajectory from which 
time-dependent NH stretching frequencies are extracted using the correlation 
between the NH fundamental transition gap and the N • ■ • N distance [22J (Sec- 
tion EH]). In Section [3721 the frequency correlation function is determined and 
lineshape theory [23] is used to obtain the IR absorption spectrum which is 
compared with the experiment results of Ref. [llj . 



2 Simulation Details 

Following the experimental setup of Ref. [11] we will investigate the WC form 
of 9-ethyl-8-phenyladenine (A) and 1-cyclohexyluracil (U) (see, Fig. [I]) sol- 
vated in CDCI3. Our model is based on the hybrid QM/MM method as imple- 
mented in the Gromacs/CPMD interface [23]. All atoms of A:U pair other than 
those in the three substituents are dealt with by CPMD [25]. The substituents 
as well as the solvent are handled by Gromacs [26] using the OPLS all-atom 
force field [27]. Since chemical bonds are cut by the QM/MM boundary, H 
atom capping is used to saturate the dangling bonds (see Fig. [TJ. 

In the simulation, one A:U pair is solvated in 100 CDCI3 molecules within a 
box with dimensions 30.0 A x 23.5 A x 23.5 A (density is 0.1 M [TT])- The 
QM part is placed in a 21.2 A x 15.9 A x 15.9 A box and the Becke exchange 
and Lee- Yang-Par correlation functional (BLYP) together with the plane wave 
basis set is used as implemented in the CPMD code. Further the Vanderbilt 
ultrasoft pseudopotential is employed with a wave function cutoff of 30 Ry. The 
MM molecular dynamics run is performed at 298 K using a time step of 2 fs 
and the Car-Parrinello dynamics is integrated with a time step of 0.121 fs and 
a fictitious electron mass of 400 a.u.. The initial configuration corresponded 
to the geometry of the optimized gas phase base pair replacing solvent atoms 
in the equilibrated simulation box. Subsequently a 7.5 ps trajectory has been 
generated, the first picosecond of which was assigned for equilibration and not 
used for data analysis. 

It is well appreciated that there is a correlation between the HB length, i.e. here 
the N • ■ ■ N distance, and the respective NH stretching frequency. Using the 
correlation curve obtained by Novak [22] supplemented by more recent data 
we can establish the time-dependent fundamental NH stretching frequency, 
u>io(t), along the trajectory and its correlation function 

C(t) = (K)(t) - (w 10 )) (wio(0) - (wi ») • (1) 
Here, (u>i ) is the mean value along the trajectory. In passing we note that 
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frequency-distance correlations also exist for O-H(D). . .0 HBs [28]. They have 
been used, e.g., in Ref. |29j to simulate lineshape functions for HOD in D 2 0. 

Adopting an adiabatic two-level model and the Condon approximation where 
all other degrees of freedom except the NH stretching frequency are treated 
classically, the IR spectrum can be expressed via the lineshape function 



Z T 

g(t) = Jdrj dr'C(r'), (2) 



as [231 



1 00 

a(u) = -Re / dte i(w - <Wl0>) *- fl(t) . (3) 
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The correlation function, Eq. (JT]), is a classical quantity and does not satisfy 
detailed balance. In principle detailed balance can be restored by applying 
certain correction factors. On the other hand, one might ask whether this mat- 
ters at all for room temperature simulations which usually are characterized 
by rapidly decaying correlation functions. This problem has been addressed 
recently by Skinner and coworker for HOD in D2O, who showed that quan- 
tum corrections have a modest influence on the IR lineshape only, but are 
important for nonlinear spectroscopies [30]. Below we will take into account 
quantum effects by adopting the multimode oscillator model [23] for fitting 
the classical correlation function by the expression 



C(t)=C'(t)+iC"(t) 

= ^2 Sj^j coth(hujj/2kBT) cos(c<j jt) + i SjU? sin(a^t) . (4) 



Here, Sj is the dimensionless Huang-Rhys factor giving the coupling strength 
between the two- level system and the mode with frequency uuj. The real part 
C'(t) is first used to fit the correlation function obtained from the molecular 
dynamics trajectory to a set of oscillators. Subsequently, C"(t) is determined 
and the complex lineshape function is used to calculate the spectrum. 
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3 Results and Discussions 



3.1 Trajectory Analysis 



The geometry is one of the most prominent characteristics of HBs and provides 
via correlation relations access to spectroscopic observables such as fundamen- 
tal IR transitions [22] or NMR chemical shifts [21] • In Fig. [2^ we show the 
N • • • N separation along the production part of the trajectory. Apparently, 
the HB length varies over a large range, that is, from 2.62 A to 3.94 A. The 
time average is 3.2 A for the N • • • N separation and 2.2 A for the H • • • N 
distance. The latter value is at the upper boundary of what has been reported 
for purine and pyrimidine crystals, that is, 1.73 A < L(H • • • N) < 2.23 A with 
a mean of 1.882 A [32]. 

Fig. [2b presents the distance between the H atom and the geometric center 
of the N • • • N connecting line. The distance is defined to be negative if the 
H atom transfers to the adenine side. First, we notice that during most of 
the simulation interval the H atom is close to the uracil Nitrogen. Hydrogen 
transfer occurs in a short period of time only, that is, in the time interval 
6.5 ps < t < 6.56 ps. Apart from this short interval there is a clear corre- 
spondence between N • ■ ■ N and N — H ■ ■ • N distance changes. This indicates 
that the HB should be close to being linear. To clarify this point further the 
deviation of the H atom away from the linear motion in the HB defined as 
the difference between the sum L(N — H) + L(N • • • H) and the bond length 
L(N---N), is drawn in panel (c) of Fig. [21 Significant deviations occur for 
times greater than t = 6.5 ps only and can be traced back to the H transfer 
event. Here, the H atom acquires additional kinetic energy after it transfers 
back to the uracil side. The reason lies in the fact that the potential energy 
curve is an asymmetric double well with the higher energy at the adenine side. 
The additional kinetic energy will lead to a substantial out of line motion un- 
til it is dissipated to the A:U pair and the solvent, as shown in the window 
around t = 7 ps of this panel. Except for the period of H transfer and relax- 
ation thereafter the difference is smaller than 0.09 A for 95% of the simulation 
time and smaller than 0.14 A for 99% of the time. The average deviation is 
0.034 A for times before the H transfer and 0.044 A for the whole simulation. 
Compared with the average NH distance (1.05 A), the deviation is small and 
the HB can indeed be considered to be linear. 

In a next step we will correlate the N • • • N separation to the NH stretching 
frequency using the experimentally established correlation of Novak [22] which 
we supplement by two additional data points for small distances (see caption 
Fig. [3]). The frequency- distance correlation is shown in Fig. [3h together with 
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a linear regression fit to the function 



f(r) = u;oo/27rcerf ^ a * r ' 



(5) 



i=0 



where Uoo is the free NH stretching frequency in gas phase (for parameters see 
figure caption). As a note in caution we emphasize that this type of correlation 
curve is of approximate nature. On the experimental side, the spectrum is 
often very broad and bears sub-band structure which makes the assignment 
difficult. On the computational side, this model assumes that the H atom is 
attached the original site, therefore, there is no H transfer at all. Of course, 
this is not the case throughout our simulation and we will have to exclude the 
corresponding configurations. Further we note that Eq. (jSj) also should not 
be used for very strong HBs, since there the transition frequencies are very 
sensitive to the details of the particular potential energy surface as well as to 
zero point energy effects (see, e.g. Ref. [33]). 

With help of the frequency- distance correlation, the transition frequency evo- 
lution along the trajectory can be obtained from the HB N • • • N length. The 
result is shown in Fig. [3b which also displays the rescaled HB length to em- 
phasize the correlation between stretching frequency and HB length. Notice 
that the correlation curve levels off at large N • • • N distances such that the fre- 
quency essentially stays constant. The average frequency along the trajectory 
in the interval 1.0-6.2 ps is 3227 cm' 1 . 

3.2 IR Lineshape 

Having calculated the time-dependent frequency u>io(t) we can proceed to de- 
termine the correlation function, Eq. ([1]). This is done by assuming ergodicity, 
that is, C(t) is obtained along the trajectory via 



where [Tj,Tj] is a certain propagation time interval. Concerning this time 
interval one has to be cautious to ensure that the time Tj is chosen such 
that the system has equilibrated from the initial configuration. The impact of 
different choices of the time interval is shown in Fig. and b. One finds that 
when Tj varies from Tj = 1.0 ps to Tj = 1.4 ps the differences of the correlation 
function are minor for t < 1.8 ps and quite small for t > 1.8 ps. This finding 
indicates that as far as the simulation of C (t) is concerned the system is indeed 




dT(uj w (t + r) - (u w ))(uj w (t) - (u w )) 



(6) 
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sufficiently equilibrated at t = 1.0 ps. Panel (b) shows that for Tj = 1.0 ps, 
the result for Tf = 6.0 ps is almost identical to that for Tf = 6.2 ps but 
quite different from that obtained using Tf = 7.2 ps. The reason lies in the 
previously mentioned fact that the H transfer which occurs at about t = 6.6 ps 
leads the system away from the uracil NH equilibrium site. In other words, H 
transfer is a singular event in the considered time interval and in order to make 
any statement concerning its influence one should have at hand a much longer 
trajectory sampling many of such H transfer events. In addition the correlation 
in Fig. [3^ has been developed for local NH stretching vibrations and becomes 
questionable if the N-H bond is replaced by a N- • -H bond. Combining the two 
panels, one can draw the conclusion that the correlation function after time 
t = 1.8 ps is sensitive to the time interval and is not reliable. Based on this 
consideration, the time range 1.0 ps < t < 6.2 ps will be used for the time 
average. 

To proceed we will fit C{t) using the multimode oscillator model, Eq. (TjJ, as 
outlined in Section [2j The frequency interval is dictated by time-step and in- 
terval of the correlation function, i.e. 80-1770 cm -1 , discretized into 20 modes 
whose frequencies and coupling constants are fitting parameters. Given the 
number of parameters, the correlation function can be fit almost perfectly in 
the interval up to 1.8 ps. The fitted coupling strength is a smooth function 
of the frequency, as shown in Fig. [4b. Specifically we find that the coupling 
between the NH vibration of the A:U pair and all other degrees of freedom 
is dominated by low-frequency modes. This includes A:U modes modulating 
the HB length (N---N distance). Preliminary normal mode calculations for 
the gas phase A:U pair show that the shoulder around 400 cm -1 could also be 
due to various modes having N ■ • • N stretching character. In passing we note 
that this shoulder is reproducible for different sets of starting frequencies in 
the fitting. 

Finally, in Fig. [4H we compare the simulated and measured [UJ IR absorption 
spectrum. The experimental spectrum peaks at 3185 cm -1 and has a full-width 
at half maximum (FWHM) of 53 cm -1 . The average transition frequency of 
our simulation has been 3227 cm -1 which upon incorporating the harmonic 
correction factor shifts to 3209 cm -1 , i.e. the difference with respect to the 
experiment is only 1%. The calculated FWHM is 39 cm" 1 which is 24% smaller 
than the experimental value. Both spectra show some asymmetry with respect 
to the high energy side. We note that without the quantum correction the 
linewidth reduces to 32 cm -1 . Besides the perhaps surprisingly good agree- 
ment between theory and experiment the analysis in Fig. [4H allows to assign 
the spectrum of modes coupled to the NH stretching vibration. One notices 
that this spectrum is cut-off at high frequencies. Therefore, contributions from 
the first overtone of the NH bending vibration are not explicitly included. It 
can be argued that possible Fermi resonances are implicitly included in the 
frequency-distance correlation. However, to establish the role of such a Fermi 
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resonance a model which includes both vibrations, e.g., by calculating respec- 
tive potential energy surfaces for snapshots along the trajectory, should be 
used. 

In summary we have presented a QM/MM approach to the calculation of 
the IR lineshape of the H-bonded NH stretching vibration of an A:U pair in 
CDCI3 solution. The time-dependence of the transition frequency was obtained 
using the N • ■ ■ N-distance-frequency correlation from Novak [22] supplemented 
by more recent data and fitted to a function which contained the free NH 
vibration as a limiting case. The obtained average N • ■ ■ N distance of 3.2 A 
corresponds to a transition frequency of 3344 cm -1 , which is slightly larger 
than the experimental value (3185 cm -1 [H]). However, using lineshape theory 
together with a multimode oscillator model an IR spectrum has been obtained 
which is quantitative as compared to the experiment within 1% for the line 
position and 24% for the FWHM. The spectral density responsible for the 
width of the spectrum is dominated by modes in the range below 200 cm -1 . 
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Fig. 1. Structure of Watson-Crick pair formed by 9-ethyl-8-phenyladenine and 1-cy- 
clohexyluracil and separation into QM and MM regions. The three H atoms bonded 
to N9 and Cg of adenine and Ni of uracil, respectively, are capping atoms in the 
QM/MM treatment. 

Fig. 2. Geometric parameters of the N — H • • • N hydrogen bond along the tra- 
jectory: (a) HB length; the horizontal line indicates the time average, (b) distance 
between the H atom and the center of the N • • • N connecting line, (c) out-of-line mo- 
tion of the H atom as measured by the difference L (N — H) + L (N • • • H) — L (N • • • N) . 
Notice that the first picosecond of the trajectory has been assigned for equilibration 
from the initial configuration. 

Fig. 3. (a) Empirical frequency-distance correlation function. The dots give the 
frequencies obtained from crystals containing intermolecular N — H • • • N HBs by 
Novak [22] which are supplemented by the points for 4-aminopyridine hemiperchlo- 
rate (2240 cm" 1 , 2.698 A) [M] and phtalazine hemiperchlorate (2370 cm" 1 , 2.656 
A) [35]. The solid line corresponds to the fitting function in Eq. ([5]) using the free NH 
stretching frequency of uracil in the gas phase {uJoo/2-kc = 3436 cm -1 |36j). The fit 
parameters for Eq. ([5]) are ao = 7.0911, a\ = —5.7941 A , and a<i = 1.2711 A . (b) 
Transition frequency o^io as a function of time. Solid line: frequencies fitted from 
the frequency-stretching correlation function, Dashed line: HB length L(N • • • N) 
shifted and rescaled such as to enable correlation with the frequency, Dotted line: 
3216 cm" 1 , the time average of the frequency along the trajectory. 

Fig. 4. Frequency fluctuation correlation function, Eq. ([6]), and IR spectrum, Eq. 
([3]). (a): impact of different choices for initial time, (b): impact of different choices 
for the end of the time interval, (c): coupling strength between the NH vibration 
and the oscillator modes as a function of their frequency (points are fitted results 
and line is guide for the eye), (d): IR lineshape of NH stretching compared with 
experimental data [11] (line for experimental result is drawn as guide for the eye 
only). 
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Fig. 2 (Yan et al.) 
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Fig. 4 (Yan et al.) 
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